# -*- coding: utf-8 -*-

import numpy as np
from datetime import datetime, timedelta

class Algorithm:
    # 定义常量
    Gsc = 118.108
    a = 0.248
    b = 0.752

    def get_algorithm(self,doy:datetime,lat:float,illumination:float):
        """
        计算光合有效辐射PAR
        :param doy: 日期
        :param lat: 纬度
        :param illumination: 光照
        :return: float PAR
        """
        # 处理日期
        doy2 = datetime(doy.year,1,1)
       
        # 计算 DOY
        days = (doy - doy2).days + 1
        
        # 处理经度
        LAT = lat * np.pi / 180
            
        # 计算公式
        F = 2 * np.pi * (days - 1) / 365
        V = (0.006918 - 0.399912 * np.cos(F) + 0.070257 * np.sin(F) -
            0.006758 * np.cos(2 * F) + 0.000907 * np.sin(2 * F) -
            0.002697 * np.cos(3 * F) + 0.00148 * np.sin(3 * F))
        
        Ws = 1 / np.cos(-np.tan(LAT) * np.tan(V))
        Ws1 = np.arccos(-np.tan(LAT) * np.tan(V))
        
        SL = (2 / 15) * Ws * (180 / np.pi)
        E0 = (1.00011 + 0.034221 * np.cos(F) + 0.00128 * np.sin(F) +
            0.000719 * np.cos(2 * F) + 0.000077 * np.sin(2 * F))
        
        H0 = (1 / np.pi) * self.Gsc * E0 * (np.cos(LAT) * np.cos(V) * np.sin(Ws) +
                                        np.sin(LAT) * np.sin(V) * Ws)
        
        HL = 0.8 * H0
        H = HL * (self.a + self.b * 0.1 * illumination / SL)
        PAR = 0.5 * H
        
        return PAR
